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We present results from a numerical study of critical gravitational collapse of axisymmetric dis- 
tributions of massless scalar field energy. We find threshold behavior that can be described by the 
spherically symmetric critical solution with axisymmetric perturbations. However, we see indica- 
tions of a growing, non-spherical mode about the spherically symmetric critical solution. The effect 
of this instability is that the small asymmetry present in what would otherwise be a spherically 
symmetric self-similar solution grows. This growth continues until a bifurcation occurs and two 
distinct regions form on the axis, each resembling the spherically symmetric self-similar solution. 
The existence of a non-spherical unstable mode is in conflict with previous perturbative results, and 
we therefore discuss whether such a mode exists in the continuum limit, or whether we are instead 
seeing a marginally stable mode that is rendered unstable by numerical approximation. 



I. INTRODUCTION 

In this paper we present results from a numerical study 
of critical collapse of the massless scalar field in axisym- 
metry. In spherical symmetry, the threshold of black hole 
formation was first systematically explored in PJ, which 
described intriguing behavior, called critical phenomena, 
in solutions approaching the threshold. This behavior 
includes power-law scaling of the mass M of black holes 
that form in the super-critical regime, 

Mcxb-p*) 7 , (1) 

where 7 is a universal constant (i.e. independent of the 
initial data) called the scaling exponent. Here, p is a pa- 
rameter describing some aspect of the initial distribution 
of scalar field energy such that for p > p* black holes 
form during evolution, while for p < p* all of the scalar 
field disperses to infinity. Thus, p* denotes the threshold 
of black hole formation for the particular family of ini- 
tial data under consideration. The solution approached 
in the limit p — > p* , called the critical solution, was also 
conjectured to be universal, in that all one-parameter (p) 
families of initial data having a threshold parameter p* 
should exhibit the same critical solution in the vicinity 
of collapse. In addition, the critical solution for the real 
scalar field is discretely self-similar, characterized by an 
echoing exponent A. Since the initial discovery reported 
in 0, critical phenomena have been observed in numer- 
ous systems — see 0, @] for recent review articles on the 
subject. Note that the particular behavior observed in 



the threshold solution depends upon the matter model 
and spacetime dimensionality. 

To date, the only non-perturbative calculation of crit- 
ical gravitational collapse away from spherical symmetry 
was carried out by Abrahams and Evans Q , who studied 
the collapse of pure gravitational waves in axisymmetry 
(note that axisymmetry is the 'minimal' symmetry one 
can impose on gravitational waves and retain the possi- 
bility of black hole formation) . In addition, the thresh- 
old of singularity formation in a non-linear sigma model 
in three dimensions was considered in [5], and found to 
exhibit features similar to critical gravitational collapse. 
These studies provide evidence that critical phenomena 
is observed beyond spherical symmetry at the respective 
thresholds of these two distinct physical systems. 

An explanation for critical phenomena, in particular 
the observed universality of the solution and departures 
from it in near-critical collapse, is offered by positing that 
the critical solution, when perturbed, has exactly one un- 
stable mode 6]. That there is only one unstable mode 
allows the threshold solution to be found in a numerical 
collapse "experiment" whereby we fine-tune a single pa- 
rameter of a generic family of initial data. Furthermore, 
the nature of the unstable mode eventually dominates 
the properties of near-critical solutions; for example, the 
scaling exponent 7 can be shown to be equal to the in- 
verse of the exponential growth factor A of the unstable 
mode. 

The purpose of the present study is to move beyond 
spherical symmetry and to explore the threshold of black 
hole formation from the collapse of axisymmetric distri- 



2 



butions of the massless scalar field. Linear perturbation 
studies of the scalar field critical solution beyond spher- 
ical symmetry were carried out by Martin-Garcia and 
Gundlach (a similar analysis has also been performed 
by Gundlach Q for the case of perfect fluid collapse). 
Their study found no additional growing modes beyond 
the one seen in spherical collapse, a result that suggests 
that we should expect to see the spherical critical solution 
emerge from our axisymmetric studies. 

Having looked at a variety of initial configurations of 
the scalar field, some deviating significantly from spher- 
ical symmetry, we find that in all cases during the early 
phases of near-critical evolution, we do see a discretely 
self-similar solution unfold that can be described as the 
spherically symmetric critical solution plus perturba- 
tions. However,in contrast to the perturbation theory 
calculations in |7j, we find some evidence for a second, 
slowly growing unstable mode, with an angular depen- 
dence described by the 1 = 1 spherical harmonic. The 
simulations suggest that this mode will eventually cause 
a near-critical solution, with some asymmetries, to "bi- 
furcate" into two distinct echoing solutions, which, in- 
dividually, would subsequently be subject to the same 
instability. In principle then, if we could fine tune to 
arbitrary precision, this bifurcate behavior would be re- 
peated indefinitely. 

The appearance of this second unstable mode is in 
conflict with the above-mentioned pcrturbative results. 
One possibility is that the non-spherical mode that ap- 
pears unstable in our simulations is in fact damped in the 
continuum limit, and only grows within the context of 
our discrete numerical approximation. Our current code 
(running on the computer systems to which we have ac- 
cess) cannot provide the accuracy needed to conclusively 
determine that the growth rate of the suspect mode is 
positive in the continuum limit, and not dominated by 
truncation error effects. 

The remainder of this paper is organized as follows. In 
Sec. |Ii|we briefly describe the relevant system of equa- 
tions, the numerical code used to solve them, and various 
properties of the solution that we will analyze. Details of 
the formalism and numerical technique can be found in 
. In Sec. IIIII wc describe several of the families of initial 
data that we have studied, and present the results from 
corresponding near-critical collapse simulations. We con- 
clude in Sec. HVI bv summarizing the results and possible 
future directions of study. 

II. PHYSICAL SYSTEM AND ANALYSIS OF 
SOLUTION PROPERTIES 

We are interested in solving the Einstein field equations 

R^u - i^Rg^u = 8ttT^, (2) 

where R^ u is the Ricci tensor, R = R p „ is the Ricci 
scalar, and we use geometric units with Newton's con- 



stant G and the speed of light c set to I. We adopt a 
massless scalar field $ as the matter source, with corre- 
sponding stress-energy tensor given by 

7^ = 2$^-^$^, (3) 

and the evolution of $ is governed by the wave equation 

□$ = $./ = 0. (4) 

Note that © differs by a factor of 2 from the convention 
of Hawking and Ellis [Hj , which amounts to rescaling $ 
by a factor of 

We restrict our attention to axisymmetric spacetimes 
without angular momentum, and choose the following 
cylindrical coordinate system, adapted to the symmetry 

ds 2 = -a 2 dt 2 + i; A [ {dp + (3 p dt) 2 

+ (dz + (3 z dt) 2 + p 2 e 2ps d4) 2 }. (5) 

The axial Killing vector is (d/d4>) p and hence all the 
metric functions a, (3 P , (3 Z ,ip and a, and scalar field <& 
depend only on p, z and t. 

We use the (2+I) + I formalism to arrive at the sys- 
tem of partial differential equations (PDEs) that we need 
to solve, which in the absence of angular momentum is 
the same set of PDEs that the ADM decomposition pro- 
vides. The Hamiltonian constraint yields an elliptic PDE 
for the conformal factor ip, and the p and z momentum 
constraints give elliptic PDEs for the p and z components 
of the shift vector, j3 p and /3 2 , respectively. We choose 
maximal slicing, in particular K a a = 0, where K a h is the 
extrinsic curvature tensor of t = const, slices; this con- 
dition gives an elliptic equation for the lapse function a. 
We convert the hyperbolic evolution equations for a and 
$ to first order form by defining "conjugate" variables f2 
and n by 

- -1KS-KS 

n = p - — (6) 

p 

and 

n = ^($. t -^$. p + r$ z ) ( 7 ) 

a 

respectively. 

We thus end up with a mixed hyperbolic-elliptic sys- 
tem of PDES for the 8 variables a, -0, P p ', (3 z ,tl,& and 
II that we approximately solve using second-order accu- 
rate finite difference (FD) techniques. The hyperbolic FD 
equations are solved using an iterative Crank-Nicholson 
scheme with adaptive mesh refinement (AMR), and the 
elliptic FD equations are solved (on the adaptive grid hi- 
erarchy) using the FAS multigrid algorithm. At t = 0, we 
freely specify a, f2, $ and II, then solve the 3 constraint 
equations and slicing condition for the remaining vari- 
ables. After t = 0, we continue to use the momentum 
constraints to solve for (3 P and f3 z , and the slicing con- 
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dition for a, but in lieu of the Hamiltonian constraint, 
we update ip using the first order in time evolution equa- 
tion that follows from the definition of the the extrin- 
sic curvature. (Thus, we employ a partially constrained 
evolution.) We add Kreiss-Oliger dissipation to the dif- 
ferenced form of the hyperbolic equations, to reduce un- 
wanted (and un-physical) high-frequency components in 
their solutions. For outer boundary conditions, we apply 
outgoing radiation (or Sommerfeld) conditions on <f>, n, a 
and Q, and appropriate asymptotic fall-off behavior for 
the remaining variables, assuming an asymptotically flat 
coordinate system. 

More details on the boundary conditions, system of 
equations and the numerical scheme (including various 
tests) can be found in 0; a detailed description of the 
AMR implementation is given in [T^ ] 



A. Analysis of Solution Properties 

In Section [illl we will quantitatively describe the near- 
critical solution for any given family of initial data by 
measuring its associated scaling exponent (7), echoing 
parameter (A), the local minima/maxima attained by 
the scalar field during each half-echo, and deviations of 
the scalar field from a spherically symmetric profile. That 
we can define such properties for all the solutions is an in- 
dication that they are similar enough that a comparison 
is meaningful. However, it is not a trivial task to compute 
some of these quantities, because we need to make sure 
that we are calculating them in a coordinate indepen- 
dent fashion. Our coordinate system is not "symmetry- 
seeking" ^jj, and the initial data is sufficiently different 
among the various families that we can expect, and in 
some cases clearly see, "gauge" differences between solu- 
tions that are apparently quite similar. 

The simplest quantity to calculate is the scaling expo- 
nent 7. We use the method proposed by Garfinkle and 
Duncan [T^j. whereby we measure the maximum value 
i? m attained by the absolute value of the Ricci scalar, \R\, 
in a set of sub-critical evolutions; 7 can then be obtained 
from the following property of near-critical solutions 

ln|i? m | w ~2-f\np + w(lnp) + const. (8) 

Here p = p* — p and w is a periodic function of its argu- 
ment with period A/ (27) that describes a small "wiggle" 
superimposed on an otherwise linear relationship. As a 
result, we can also use (JSJl to obtain an estimate for A. 
We note that the effectiveness of our use of JHJ) to com- 
pute 7 and A is predicated on the degree to which our 
computed near-critical solutions are well approximated 
by a discretely self-similar solution with a single unsta- 
ble mode. In addition, although © provides the only 
method we use to estimate 7, we also measure A using a 
more direct procedure outlined below. 

The direct comparison of results from our axisymmet- 
ric code to those from a spherically symmetric computa- 



tion presents more of a challenge. In order to compare 
"local self-similar solutions" — portions of the computed 
spacetime that appear to be approximately self-similar 
about some center of symmetry — we need an invariant 
way of slicing the spacetime in the region of interest. 
To accomplish this, we use a sequence of outgoing null 
hypersurfaces, starting from the local center of symme- 
try (p, z) = (0, zq), to generate the common slices along 
which we compare <J>. To construct each such null hyper- 
surface, we evolve a family of null geodesies, with affine 
parameter x and initial tangent vectors equally spaced 
in 9 = tan _1 (/9/(z — zq)), outward from (p,z) — (0,zq). 
The geodesies are synchronized by setting dx/dr = 1 at 
the start of integration, where r is the proper time mea- 
sured by a timelike observer that is stationary relative 
to the center of symmetry, r is the time of relevance to 
critical collapse, for in coordinates ln(r* — r), where t* 
is the accumulation point of the critical solution (i.e. the 
central proper time of the central singularity formed by 
the cascade of the critical solution down to infinitesimally 
small scales) , the central value of the scalar field is a pe- 
riodic function of ln(r* — r), with period A. Estimation 
of the period of the profile of <E> along the local center of 
symmetry, with respect to r thus gives us the alternate 
method for computing A. 

During a simulation, we integrate function of 

t for each null geodesic labeled by 80 = 9(x = 0), and 
record $(2;, 9q) (we typically use 50 geodesies per slice, 
linearly spaced in 9q). If two solutions from different fam- 
ilies of initial data do locally tend to the same discretely 
self-similar solution, then $(x,0o) (synchronized so that 
the null integration is started at the same time within 
the periodic oscillation) will tend to the same function, 
regardless of differences in the (p, z, t) coordinate systems 
between the two solutions. 

As a final comment in regards to our analysis, we 
note that to calculate r we integrate a central timelike 
geodesic, and measure proper time along it. We can do 
this for families of initial data that are symmetric about 
2 = 0, for then we know that the center of symmetry will, 
at least initially, be at (0,0). An interesting aspect of the 
numerical solution is that truncation error effects cause 
a small drift to occur in the z location of the local cen- 
ter of symmetry, during a near-critical evolution. This 
drift is quite small (and does appear to converge away 
with increasing resolution), typically being less than 1 
part in 10 6 of the size of the computational domain. 
However, because of the exponentially decreasing length 
scales that arise in a critical collapse, this is a huge drift 
relative to the size of the local self-similar region at late 
times. Hence, if we simply measured central proper time 
at (p,z) = (0,0), and correspondingly integrated null 
geodesies from this location, we would entirely miss the 
relevant part of the solution. Fortunately, the timclikc 
observer initially placed at (0, 0) experiences an identical 
drift, and so we can use its location and proper time to 
do the desired measurements. For initial data that is not 
plane-symmetric (the only such family described in the 
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next section is the "anti-symmetric" example) we have 
not yet been able to devise a method to accurately track 
the local center of symmetry for long periods of time, 
and hence have not been able to calculate t for these 
families. However, at least at key moments during the 
evolution, we are able to accurately determine the center 
of symmetry (by looking at local minima or maxima of 
<fr, for instance), to use as the starting point for the null 
integration. 



III. RESULTS 

Here we present results from the critical collapse of 
several families of initial data. These families consist of 
a time-symmetric series of prolate spheroids, with ellip- 
ticity e (defined by looking at surfaces of constant <&): 

$(0,p,z) = Ae-( p2+{1 ~ e2)z2 ), 

n(o, P ,z) = o, (9) 

and an initially ingoing distribution in $ that is anti- 
symmetric about z = 0, i.e. $(<, p, z) — — ■!>(£, p, —z): 

H0,p,z) = Aze-i^^-* ) 2 , 

n(o, P ,z) = -$(o,p,z). (io) 

In all cases we set a(0,p,z) = and O(0, p,z) — 0, and 
vary the amplitude A when searching for the threshold 
solution. We show results for six families of ©, with 
e 2 = 0,1/3,1/2,2/3,3/4 and 5/6. In R is a pa- 

rameter describing how far the initial pulse of matter 
begins from the origin; we have chosen i?o = 3. In all 
cases presented here the outer boundary of the compu- 
tational domain is at p — \z\ — 10, though in other sim- 
ulations we have varied its position to make sure that 
the above choice does not significantly impact the results. 
The base level in the adaptive hierarchy used a resolution 
of 65 x 129 points, and up to 28 additional 2:l-refined lev- 
els were used in the most nearly critical case. Our AMR 
implementation is based on the algorithm of Berger and 
Oliger 15], wherein regridding is determined through es- 
timates of the local truncation error (solution error) in 
the computed solution. The key control parameter that 
determines placement of refinements is the truncation er- 
ror threshold, r m : mesh refinements are introduced in an 
attempt to keep the magnitude of the local truncation er- 
ror estimate < r m throughout the solution domain. For 
each family of initial data studied, we generally tuned 
to threshold using three different values of r m , namely 
r m() , T m0 /2 and r m0 /4. Most of the data presented here 
are from r m = r m o/4 runs (i.e. finest effective resolu- 
tion), with the results computed using the less stringent 
values of r m then being used to give some estimate of 
how close to the continuum solution we may be (though 
convergence testing with an adaptive code is not trivial, 
particularly in the critical limit). 



There are of course, infinitely many different 
parametrized families that we could have considered — 
those used to generate the results discussed here were 
chosen for the following specific reasons. First, the anti- 
symmetric configuration l|10|) provides a more drastic de- 
parture from spherical symmetry than any family of data 
that can smoothly be deformed into a spherical distribu- 
tion (such as © by letting e — > 0). In this regard we 
note that one of the characteristic features of spherical 
scalar field critical collapse is that the "central" value of 
$ oscillates between specific extremal values ±$o; clearly 
the anti-symmetric property of l|10l) allows no such oscil- 
lation. One might therefore expect that evolutions with 
this type of initial data might produce a qualitatively 
different critical solution than the spherically symmet- 
ric one. However, as shown in Figure ^ a t threshold 
two spherical-like echoing solutions develop off-center at 
z = ±z c (t). 

Second, we include the prolate family because the 
initial amplitude of the putative second unstable mode 
seems to be closely related to the prolateness of the ini- 
tial distribution (rather, for instance, than asymmetries 
in «3> within an imploding spherical shell). Therefore, the 
parameter e in @ allows us to demonstrate the effect 
of adding more (larger e) or less of the unstable mode. 
The axisymmetric instability, once it has grown beyond a 
certain amplitude, causes a near-spherical threshold solu- 
tion to "bifurcate" into two echoing solutions, separated 
by some distance along the axis. As an example, Figure 
121 shows several time- instants from the near-critical evo- 
lution of initial data with e 2 = 3/4, transformed to loga- 
rithmic coordinates in space to better illustrate the self- 
similar nature of the initial critical behavior. Note that 
this bifurcation is qualitatively different from the two 
echoing solutions observed in anti-symmetric collapse — 
there, by construction, no self-similar behavior is seen 
about z — 0, and there are two (out of phase) echoing 
solutions from the beginning. Furthermore, the initial 
separation of the two (in phase) echoing solutions arising 
from a bifurcation is related to the smallest length scale 
that developed in the single, origin-centered echoer prior 
to the bifurcation. In contrast, the separation of the two 
anti-symmetric echoing solutions is related to a length 
scale in the initial data. Moreover, if there really is a 
second unstable mode, then each of the anti-symmetric 
echoers should also be subject to that instability and 
eventually bifurcate, and we do see some evidence for 
this. 

With double precision arithmetic, we are able to tune 
the initial amplitude of a given family to within a part 
in 10 15 of threshold 16], corresponding to about 3 full 
echoes of the spherically symmetric critical solution. The 
growth of the instability is sufficiently small that after 
3 echoes we do not yet see a bifurcation for e 2 < 2/3; 
for e 2 = 3/4 and e 2 — 5/6 we see a bifurcation after 
approximately 2 and 11/2 echoes respectively. 

Table [I] summarizes measurements made of the criti- 
cal parameters — namely 7, A and the amplitude of each 
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t = 2.76 



t = 5.504 



FIG. 1: Several frames of &(p, z, t) from the evolution of near- 
critical, anti-symmetric initial data HUH . The figures span the 
first several half-echoes of the local self-similar solutions, and 
the particular times shown correspond to when the scalar field 
reaches a local minima/maxima. The height of each surface 
represents the magnitude of and the coordinate domain of 
each figure is [0..2.5, —2.5. .2.5] in [p,z] (the axis p = is the 
nearest edge of each plot, and positive to negative z runs from 
left to right). 




FIG. 2: Several frames of 3>(r, 9, t) from the evolution of near- 
critical, e 2 =3/4 prolate initial data @. Here we have trans- 
formed to coordinates r = ln(^/ p 2 + z 2 + eo) — lneo (with 
eo = 2xl0 -4 ) and tan 9 — p/z, to give a better view of the 
initial self-similar nature of the solution, [f, 9] ranges from 
[0.. ~ 10.8, 0..7r], with the axis p = being the nearest edge 
in each figure. The height of each surface represents the mag- 
nitude of $. The times shown correspond to the times when 
$ reaches a local minima/maxima, demonstrating the bifur- 
cation that occurs after about 2 self-similar echoes of the field. 
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e 2 


7 


(A) 1 


(A) 2 


(l*c|) 





0.382 ± 2% 


3.44 ± 1% 


3.49 ± 3% 


0.431 ± 2% 


1/3 


0.380 ± 2% 


3.41 ± 1% 


3.43 ± 3% 


0.431 ± 2% 


1/2 


0.375 ± 3% 


3.37 ± 1% 


3.39 ± 4% 


0.430 ± 2% 


2/3 


0.346 ± 3% 


3.13 ± 1% 


3.08 ± 4% 


0.419 ± 3% 


3/4(a) 


0.313 ±4% 


2.87 ±4% 


3.03 ± 5% 


0.396 ± 6% 


3/4(6) 


0.40 ± 10% 




w 3 


0.40 ± 8% 


5/6(a) 


0.28 ± 10% 


« 2 


« 1 


0.36 ± 7% 


5/6(6) 


0.41 ± 10% 




w 3 


0.36 ± 10% 


AS 


0.383 ± 2% 




3.49 ± 3% 


0.434 ± 3% 




-1ii(t-t*) 



TABLE I: Critical parameters of the prolate families © and 
the anti-symmetric family (AS) HUH . 7 is obtained from a 
least-squares fit to the data shown in Figure 

H (A) 1 is the 

average value of A measured between adjacent extremes in 
"1> C as shown in Figure |S] (only using data from intermedi- 
ate times), (A) 2 is the average value of A inferred from the 
periodic oscillations in Figure 01 and (|$ c |) is the average ab- 
solute value of the extremes of $ c in Figure (again using 
data from intermediate times) . For the two prolate cases that 
bifurcate — e 2 = 3/4 and e 2 = 5/6 — we list estimates of these 
parameters (where possible) before (a) and after (b) the bi- 
furcation. For the anti-symmetric case, we do not have data 
for $ c versus central proper time; (|3> e |) in that case is calcu- 
lated as half the average difference between subsequent local 
extremes in <!>(p = 0,z,t) about one of the local self-similar 
solutions. See the text for a discussion on how the estimated 
uncertainties were calculated. 



echo in <E> — from the r m o / 4 simulations for each family of 
initial data (except for the e 2 = 5/6 case, where the in- 
creasing computational demands, resulting from larger, 
more elongated grids that are produced in the hierarchy 
for higher values of e, prevented us from computing with 
anything but r m — T m o). For the two simulations with 
the largest values of e, we list parameters obtained be- 
fore and after the bifurcation, where possible. As with 
the anti-symmetric case, our method of geodesic inte- 
gration cannot track moving centers, and so we cannot 
provide a direct estimate of A after a bifurcation. Also, 
for the e 2 = 5/6 case, we do not see a very distinctive pe- 
riodic oscillation in the lni? m vs. hip plot, and thus can 
only provide a rough guess for A from that information. 
Most of the data in this table was gathered from Figures 
0]and[5l which show lni? m vs. hip (with the linear re- 
lationship from the spherical family subtracted to better 
differentiate the plots) and $ c , the central value of <E>, vs. 
logarithmic central proper time for the prolate families 
prior to bifurcation, respectively. Also, Figure [3] shows 
the same type data displayed in Figure but for the 
case e = 0, and with the addition of an overlay of data 
obtained with a spherically symmetric ID code lj. The 
good agreement between the results from the axisym- 
metric and spherical computations provides a measure of 
confidence in the correctness and level of accuracy of our 
2D code. 



FIG. 3: Comparison between results from our 2D axisymmet- 
ric code (e = initial data) and a ID spherically symmetric 
code [lj. Plotted is the central value of $ as a function of 
logarithmic central proper time, for the most nearly critical 
simulations obtained in either case. 



The quoted uncertainty of a given value in Table ^ 
was calculated as the sum of the estimated truncation 
error from a convergence calculation using the different 
r m runs, and the standard deviation from the relevant 
averaging/fitting operation (except for e 2 = 5/6, where 
we could not estimate the truncation error as we only 
have data from a single value of r m ). However, in a sense 
these uncertainties are "optimistic," for we have not ac- 
counted for possible systematic errors. Chief among these 
(in particular away from spherical symmetry) are the as- 
sumptions of discrete self-similarity, which was used to 
define t* in Fig[5j and the assumption that the linear 
and periodic parts of Fig^Jare directly related to 7 and 
A respectively. For several of the simulations we have 
checked that the following numerical parameters are not 
significant sources of systematic error: outer boundary 
location, Dirichlet vs. Neumann conditions on cv,/3 p and 
(3 Z at the outer boundary, and free vs. constrained evo- 
lution for ip. 

Note that our coordinate system is not adapted to 
spherical symmetry, and during an evolution of e = ini- 
tial data, spherical symmetry is only preserved to within 
an amount proportional to the truncation error, and 
so will eventually exhibit the apparent second growing 
mode. In a certain sense this is a desirable feature, for at 
late times during an e = evolution this mode is the only 
one (apart from the unstable spherical mode) that should 
be visible perturbing the spherical solution; an e > evo- 
lution exhibits a host of additional, decaying asymmetric 
modes that prevent us from easily measuring the proper- 
ties of the non-spherical growing mode. To this end, in 
Figure El we show plots of the maximum absolute value 
of the I — 2 (to = 0) spherical harmonic component of 
denoted $^2, in near-critical e = collapse, as mea- 
sured along outgoing null slices of the spacetime (in other 
words, we decompose <f>(x,0o), constructed as described 
in Sec. IHJ into its spectral coefficients for each x — the 
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I = 2 component is <&t2-)- We show results from Simula- error estimate r m , demonstrating the expected behavior 
tions with 3 different values of the maximum truncation that $£2 — ► in the limit r m — > 0. 



We now argue that Fig. also gives some evidence 
that the instability we do see in the numerical solution 
may be an actual feature of the continuum solution, and 
not a truncation-error-driven phenomenon. 

Assume that the numerical solution has a well-behaved 
Richardson expansion. Then we expect any well-defined 
continuum property of the solution, such as the growth 
rate, Aj, of a perturbative mode to have a similar expan- 
sion: 

A, = Ai + hf{x) + 0(h 2 ), (11) 

where Ai is the numerically measured growth rate, h is 
the discretization scale (we have assumed a first order 
accurate discretization), and /(x) is some function of the 
continuum solution variables x. Of course, in an adap- 
tive scheme there is no single scale h; however, individual 
grids within the hierarchy do admit Richardson expan- 
sions, and hence we can loosely think of 111|) holding 
over the hierarchy with some effective h that would be 
related to the maximum truncation error estimate r m . 
In 0, the real part of the largest eigenvalue of any non- 
spherical mode perturbing the critical solution was found 
to be A2 ~ —0.02; i.e. a decaying mode, and the corre- 
sponding eigenfunction had the angular dependence of 
the i — 2,m = spherical harmonic Y^. This magni- 
tude of decay is about 100 times smaller than the growth 
rate of the dominant spherically symmetric mode, that 
has A = I/7 w 2.7. Thus, looking at Qllfl. it is cer- 
tainly plausible that in a numerical scheme, even if h were 
small enough to reasonably accurately model the domi- 
nant feature of a solution (as we evidently are from the 
comparison in Figure it might still be large enough to 
significantly affect sub-dominant features of the solution, 
such as a small Ai in i|ll|). 

We should then be able to see a significant effect when 
changing h; however, in FigureEl even though the initial 
amplitude of the asymmetry decreases as r m decreases 
(as expected), the apparent growth rate that we obtain, 
namely A2 ~ 0.1 — 0.4, does not noticeably change within 
the relatively large uncertainty of the measurement |l7j . 
On the other hand, we may still be too far from the 
convergent regime to measure A2 (so that higher order 
terms in (|llf> are still important). Note that we also can- 
not conclusively say that the growing mode we see has a 
pure I = 2 angular dependence, but it appears that the 
I = 2 mode is at least an order of magnitude larger than 
any of the other asymmetric modes we find in the spectral 
decomposition. However, it must be noted that for com- 
putations with any of the three values of r m adopted we 
use 50 points in 9 along which we integrate null curves, 
and so do not have good accuracy for determining the 



higher I modes. 

In Figure below we show a plot of the growth rate of 
§12 from e 2 = 2/3 near critical solutions. This value of 
e is the smallest, non-zero value considered that clearly 
shows growth of the asymmetry during the roughly three 
self-similar echoes of evolution; in the e 2 = 1/3 and 
e 2 = 1/2 cases, early-time evolution of the i = 2 spec- 
tral component is dominated by decaying modes. For the 
e 2 = 2/3 data in FigureQwe apparently are converging to 
the growth shown; i.e. as was the case for the spherically 
symmetric initial data, the growth does not appear to 
be truncation error dominated. Estimates of the growth 
rate from the simulation with the smallest value of r m 
gives A2 Rj 0.05 — 0.15. However, this is quite a rough 
estimate as we cannot disentangle the supposed growing 
mode from the full spectrum oil — 1 modes contributing 
to the plot shown in Figured 

Although we appear to be converging to a growth of 
the asymmetry in the e 2 = 2/3 case, and to a bifurcation 
for the e 2 = 3/4 case, this does not necessarily prove that 
the spherically symmetric critical solution has a second 
unstable mode. These families are sufficiently aspheri- 
cal that one can imagine that the bifurcation is due to 
some artifact of the initial data — in particular a "focus- 
ing" effect, as the wave front of an imploding, prolate 
distribution of the scalar field will tend to focus to two 
locations on the axis, above and below the origin. If this 
is the case though, it is rather surprising that we see self- 
similar collapse occur about a single center 'prior to the 
bifurcation. 

Finally, in Figure 00 we show comparisons of the "ra- 
dial" profile of the £ = 2 spherical harmonic component 
of <f>, measured along an outgoing null geodesic at (ap- 
proximately) the same time within a self-similar echo|l8fl. 
for the e 2 = 0, e 2 = 2/3 and anti-symmetric near-critical 
solutions (with r TO = r m o/4). In the plot, the overall 
amplitude and affine distance along each null curve is 
rescaled so that the maximum amplitude is one, and oc- 
curs at one in rescaled affine time x. That the curves 
from these representative families do approximately agree 
provides additional evidence that we are seeing a unique, 
asymmetric unstable mode. 



IV. CONCLUSION 

We have presented results from a first study of scalar 
field critical collapse in axisymmetry in the fully non- 
linear regime. We find that critical phenomena is ob- 
served at the threshold of gravitational collapse of several 
families of asymmetric initial data. The critical solution 
that unfolds at threshold can (locally) be described as 
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ln(p) 

FIG. 4: Plots of the logarithm of the maximum absolute value \R m \ attained by the Ricci scalar, during sub-critical evolution, 
vs. the logarithm of the distance in parameter space p = p* — p from the (estimated) critical parameter p* , for all families of 
initial data considered (using the lowest r m data). To avoid clutter in the figure, we have placed the data from each family on 
one of two identical panels. To facilitate comparison, the line —270 hip has been subtracted from each curve, with 70 = 0.382 
(the estimated value from the e = family). Also, the intercept of each curve has been set to p = 0. The estimated linear 
relationships used to calculate 7 in Table Q are also shown in the figure. 



the spherically symmetric critical solution found in 
with asymmetric perturbations. However, in contrast to 
the results of |7|, we find some evidence that a single 
I = 2 spherical harmonic perturbation does not decay 



with time; rather it grows at a rate of roughly 1/10 the 
magnitude of the dominant, spherically symmetric un- 
stable mode. The nature of this second unstable mode 
is such that it causes a self-similar threshold solution, 
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FIG. 5: The central value of $ as a function of logarithmic central proper time, from the most nearly critical simulations of 
the prolate collapse families. The horizontal axis for each family was shifted so that the first maximum of each curve occurs 
at — ln(r — r*) = 0. r* is the accumulation point of each family (prior to the bifurcation for the e 2 = 3/4 and e 2 = 5/6 cases), 
calculated by assuming that the intermediate time behavior is discretely self-similar, and then finding the r* for each case that 
minimizes the variance in A computed between pairs of adjacent minima/maxima in $ c . To aid in comparison, dashed vertical 
lines have been drawn at intervals of 3.44 in r, which is the estimated spherical echoing exponent. 



with some asymmetry in it, to eventually bifurcate into 
two local, self similar solutions that again resemble the 
spherical threshold spacetime. If this second instability 
is indeed a property of the spherically symmetric critical 



solution, then presumably one (or both if the initial data 
has reflection symmetry) of the new self-similar solutions 
would bifurcate again, and so on, resulting in an infi- 
nite, "random walk" of bifurcations on ever decreasing 



10 




n i I i r 



CD 
N 

c6 

a 

O 

2: 



e s -0 
e 2 =2/3 

Anti- symmetric 




FIG. 6: The maximum absolute value of the £ — 2 (m = 0) 
spherical harmonic component of <E>, $12, in near-critical e = 
collapse, measured along outgoing null slices of the spacetime, 
from simulations with 3 different values of the maximum trun- 
cation error estimate r m . The graph indicates similar growth 
rates independent of the effective resolution. In order to bet- 
ter show the similarity of the growth rates, the data for the 
higher truncation error thresholds have been shifted along the 
horizontal axis which labels logarithmic central proper time 
(about —3.4 for the to/2 case and —5.1 for the to/4 data). 
Thus, at similar un-shifted times, the graph indicates that 
the amplitude of decreases with r m , as expected. 
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FIG. 7: The maximum absolute value of the t = 2 (m = 0) 
spherical harmonic component of $, <J>«, in near-critical 
e 2 = 2/3 collapse, measured along outgoing null slices of 
the spacetime, from simulations with 3 different values of the 
maximum truncation error estimate r m . The horizontal axis 
labels the logarithmic central proper time when the given out- 
going null surface intersects the origin. Note the different ver- 
tical scales when comparing this plot to the similar one for 
e = in Figure [S] (and we have not shifted the data here). 



FIG. 8: The normalized I = 2 (m = 0) spherical harmonic 
component of 3>, measured along an outgoing null geodesic 
starting from similar times within a (chosen) self-similar os- 
cillation of the e 2 = 0, e 2 = 2/3 and anti-symmetric near 
critical solutions (from r m = r m .o/4 simulations). To facili- 
tate comparison, we have rescaled the amplitude of each curve 
so that the maximum is 1, and rescaled the affine parameter 
(labeled x) along the null curves so that the first maxima of 
each is at x = 1. 



timately affects the threshold solution, is beyond the 
capabilities of our current code. First, we are using 
double precision arithmetic, and this prohibits us from 
tuning closer than 1 part in 10 15 of the threshold. Be- 
cause the echoing exponent of the spherically symmetric 
solution is (relatively speaking) so large, 1 part in 10 15 
can only give us about three, complete self-similar 
echoes. This is far from ideal when trying to estimate 
the growth (or decay) rate of a mode that may have 
an e-folding time on the order of 10-20 echoes. Second, 
we would like to achieve higher accuracy than we have 
been able to attain so far. To do this with the current 
code will require that we parallelize it because we have 
already reached the practical limits (in terms of memory 
usage and runtime) imposed by the hardware to which 
we have access. Alternatively, one could write a code 
adapted to the spherical critical solution (for instance 
using spherical polar coordinates with a logarithmic 
radial coordinate). This would allow one to obtain 
greater resolution, with a given amount of resources, 
than what we can achieve with our more general purpose 
cylindrical coordinates. Of course, spherical polar 
coordinates would not be well suited to following a 
solution beyond a bifurcation, but it should be adequate 
to study the growth or decay of perturbations. 



scales. Thus the second instability would not completely 
destroy the universal nature of generic (axisymmetric) 
critical collapse, but rather would alter it in an intrigu- 
ing, family dependent manner. 

To conclusively answer 1) whether there is a second 
unstable mode, and if so 2) how the bifurcation ul- 
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